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Abstract 

Robust and accurate schemes are designed to simulate the couphng between sub- 
surface and overland flows. The coupling conditions at the interface enforce the 
continuity of both the normal flux and the pressure. Richards' equation governing 
the subsurface flow is discretized using a Backward Differentiation Formula and a 
symmetric interior penalty Discontinuous Galerkin method. The kinematic wave 
equation governing the overland flow is discretized using a Godunov scheme. Both 
schemes individually are mass conservative and can be used within single-step or 
multi-step coupling algorithms that ensure overall mass conservation owing to a 
specific design of the interface fluxes in the multi-step case. Numerical results are 
presented to illustrate the performances of the proposed algorithms. 
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1 Introduction 



The interactions of subsurface and overland flows are an important ingredi- 
ent for a comprehensive understanding of hydrology processes. While there is 
an extensive bulk of literature devoted to the numerical study of water flows 
in single-phase and variably saturated porous media, the issue of coupling 
such flows with surface flows generated by rivers, tides or floods has received 
less attention. One of the most popular ways to couple Darcy and Stokes 
flows is through the well-known Beavers- Joseph-Saffman condition [8,29,18]. 



This condition was used for instance in [12,24] in the mathematical and nu- 
merical study of the coupling of Darcy flow with a three-dimensional non- 
hydrostatic shallow-water model. Another approach used in numerical hydrol- 
ogy (see among others [31]) considers discontinuous pressures at the interface 
and evaluates an interface flux as the pressure difference, modulated by a mul- 
tiplicative exchange coefficient depending on the soil relative permeability. A 
third approach consists of assuming both normal flux and pressure continuity. 
This means that the hydraulic head of the subsurface flow matches the depth 
of the overland flow at the interface, while the normal ground flow velocity is 
used as a source term in the governing equation of the overland flow. Exam- 
ples of studies based on this approach include coupling one-dimensional surface 
flow with vertical soil columns [30], coupling the two-dimensional Richards' 
equation with a one-dimensional kinematic or diffusive wave approximation for 
the overland flow [22,7], and coupling the two-dimensional Darcy's equation 
with one-dimensional shallow- water equations [11] . 

In the present work, we assume that the subsurface flow occurs in a variably 
saturated porous medium and that this flow can be described by Richards' 
equation, entailing in particular that there are no trapped air pockets in the 
soil; otherwise more general multi- phase models should be used [3]. Further- 
more the kinematic wave approximation is used to describe the overland flow. 
This choice is solely made for ease of exposition and more general shallow 
water models can also be used. Concerning the coupling conditions, we adopt 
the third approach described above, namely enforcing the continuity of both 
normal flux and pressure at the interface. These coupling conditions are gen- 
erally valid when the overland flow is mainly produced by exflltration from 
the soil, so that normal flux and pressure equilibrium can be expected to hold 
at all times. A different situation, which falls beyond the present scope, is that 
of a runon surface wave rapidly propagating over a dry soil. 

Many methods can be employed to discretize in space Richards' equation, 
namely flnite differences [32,9], flnite volumes (FV) [23], flnite elements (FE) 
[9,19] or mixed flnite elements (MFE) [21,6]. These methods are generally com- 
bined with an implicit Euler time scheme. An alternative approach for space 
discretization is to use a discontinuous Galerkin (DG) method. Advantages 
offered by DG methods include local (elementwise) conservation (as FV and 
MFE), high-order accuracy (as FE and MFE) and ffexibility in the use of non- 
matching meshes (as FV), in particular within multi-physics and multi-domain 
approaches. Various forms of DG methods can be used for Richards' equation 
and more generally for two-phase flows in porous media. Examples include the 
so-called Local Discontinuous Galerkin method [15,4] and the non-symmetric 
or the symmetric interior penalty DG method [20,2,5]. In the present work, we 
choose the symmetric interior penalty DG method (in short SIPG), because it 
preserves the natural symmetry in the discrete diffusion operator. Regarding 
time discretization, the common approach when working with DG methods is 
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to employ Runge-Kutta (RK) explicit schemes [10] or diagonally implicit ones 
[2]. Here, we propose instead to use a backward differentiation formula (BDF). 
We think that this approach offers several advantages, such as high-order ac- 
curacy in the time discretization, circumventing the CFL condition which can 
be very restrictive for explicit schemes when diffusion processes are dominant, 
and in general higher computational efficiency than implicit RK schemes for 
problems where the nonlinear solver is expensive. Typically, if piecewise poly- 
nomials of degree p are used in the DG method, a BDF of order (p -|- 1) can 
be employed. 



The main objective of this work is to design robust and accurate schemes 
for coupling subsurface and overland flows. While Richards' equation is dis- 
cretized by a BDF-SIPG method, the overland flow governing equation is 
discretized by a Godunov scheme and advanced in time with a different time 
step if the overland flow time scale is quite different from the subsurface flow 
time scale. Two important issues are addressed in the design of our coupling 
algorithms: 1) satisfy as accurately as possible the coupling conditions which 
impose certain specific inequality and equality constraints on the pressures 
and normal fiuxes, similarly to the boundary conditions encountered in Sig- 
norini problems, and 2) ensure overall mass conservation for the whole system 
consisting of subsurface and overland flow. This point deserves some particu- 
lar attention. Indeed, although mass conservative schemes are used for both 
subsurface and overland flows, the interface flux must be chosen appropriately 
when working with multi-step methods such as BDFs. For simplicity, we will 
discuss in detail only the design of the interface flux for the second-order BDF. 
Finally, although the material will be presented in a 2D/1D setting (that is, a 
two-dimensional subsurface flow coupled to a one-dimensional overland flow), 
the results extend naturally to the 3D/2D setting. In particular, the wet part 
of the interface is not tracked directly, but is determined at each time step by 
a cell-oriented procedure within an iterative loop that solves consecutively the 
overland and subsurface flow governing equations. 



This paper is organized as follows. In Section 2, we present the physical prob- 
lem. In Section 3, we describe the time and space discretization of the model 
problem and design the coupling algorithms for both flrst-order and second- 
order BDFs. Finally, in Section 4, we present numerical results assessing the 
performance of the proposed algorithms on three test cases. 
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2 Model problem 



2.1 The setting 



Let f2 C denote the bounded subsurface flow domain with outward normal 
unit vector tiq. The boundary of Q is divided into three parts (see Figure 1): 
I is the upper part of the boundary where overland flow can occur, W are 
lateral walls and B represents the lower part of the boundary. At any time t, 
the set I is divided into "wet" and "dry" parts X*^'* U2^'*, with 



r^'* = {x e X; h{x, t) > 0}, T^'* {x e X; h{x, t) = 0}, 



•d,t def 



(1) 



where h is the depth of the overland flow. Observe that the above partition of 
X is time-dependent. 
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Figure 1. Schematic of the computational domain with basic notation. 



2.2 Subsurface flow 



The soil is modeled as a non-deformable porous medium in which the pores 
can contain both water and air (unsaturated zone) or only water (saturated 
zone). We assume that water is incompressible and that air pressure does not 
affect the flow. The water conservation equation takes the form 

dt[em + v-v{ij) = f, (2) 

where dt denotes partial time-derivative, i/j is the hydraulic head (m), 6{ip) the 
volumetric water content (dimensionless), v{iIj) the flow velocity (ms^^), and 
/ a volumetric water source or sink (s~^). In the sequel, we assume that there 
are no volumetric sources or sinks, so that f — 0. The flow velocity depends 
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on the hydraulic head through the generahzed Darcy law 

v{^l^)^-K{,p)V{^P + z), (3) 

where K[i/j) is the (possibly tensor- valued) hydraulic conductivity {ms~^) and 
z the vertical coordinate (m). Substituting (3) into (2), Richards' equation is 
obtained in the form [28] 

dtm)] - V • + z)) = 0. (4) 

Given at each time t G [0, T], where T is the total simulation time, the partition 
{T^'*,X'^'*} of X, a Dirichlet datum defined on X™'* and a Neumann datum 
co-t, defined on T*^'*, the subsurface flow is governed by 

' dt[9{jP)]+V-v{ij)^0 
v(jjj) = -K(^jJ)V(^|J + z) 

^(■,0)=^o 

< 

v{iIj) ■na = VN 

where uat is the possibly time-dependent normal flow velocity prescribed on 
W UB and -0° the initial condition. Thus, Richards' equation is a nonlinear 
parabolic equation which degenerates into a nonlinear diffusion equation in the 
saturated zone where 9 and K are constant. Examples for the two constitutive 
laws ip I— s> Oijp) and ip ^ K{ip), which are necessary to close the subsurface 
flow model, are specifled in §4. 

2.3 Overland flow 

Water surface flows are often modeled by a simplified form of the free boundary 
Navier-Stokes equations. Assuming hydrostatic pressure, negligible vertical ve- 
locity gradients and mild variations of the free surface leads to the well-known 
shallow-water equations; see, e.g. [16] for a derivation of these equations. Ne- 
glecting turbulence effects, the equations expressing the conservation of mass 
and momentum reduce to 

dth + d^q = {v{^) - Vr) • nn, (6) 
^gh{S-J), (7) 



in Q X [0,r], 
in Q X [0,r], 
in Q, 

on (WUi3) X [0,T], 
on eX'i'*}, 
on {{x,t),x el^'^}, 



(5) 



dtq + 5x 



h 2 
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where q is the discharge (m^s~^), v{iIj) ■ the source or sink term (ms~^) 
resulting from mass transfer between subsurface and overland flows, Vj. the 
possibly time-dependent prescribed rainfall intensity (ms~^), g the gravity 
acceleration {ms~^), S the possibly space-dependent bottom slope (dimen- 
sionless) and finally J (dimensionless) results from friction effects. Note that 
the mass transfer term v{tp) ■ uq in the mass conservation equation (6) in- 
volves the subsurface flow velocity resulting from (5); infiltration occurs if 
v{ip) ■ fiQ < whereas exfiltration occurs if v{ip) ■ na > 0. The Manning- 
Strickler uniform flow formula is chosen to link J and q and assuming the flux 
to be uni-directional from left to right so that g > 0, this yields 

q = O^/^ JV2, (8) 
where K, is the Strickler coefficient of roughness (m^/^s~^). 

A common assumption is to neglect inertia and potential energy effects in (7), 
so that momentum balance is governed by the equilibrium between slope and 
friction, that is 

S = J. (9) 

Substituting (9) into (8) yields 

q = ip{h,S)''^'lCh'/'S'/^. (10) 

Finally, using (10) in (6) leads to the so-called kinematic wave approximation 
[26] 

dth + d^(p{h,S) ^ {viip) -Vr) -uq. (11) 

This scalar conservation law is strictly hyperbolic wherever h > 0. In the 
present case, waves travel rightwards and an upstream boundary condition in 
A (see Figure 1) must be set. Let /i° be the initial condition and let /ia be the 
upstream boundary condition on the surface water depth prescribed at point 
A. Then, the overland flow is governed by 

dth + d^(p{h, S) = {v{iIj) - v^) - no, on X x [0, T], 
< /i(-,0) = /iO onX, (12) 

/i(A, •)=hA atAx[0,T]. 

2.4 Admissible set 

We refer to the quadruplet {X'^'*,X'^'*,a;^,a;^} as the coupling variables. The 
model problem considered hereafter for coupling subsurface and overland flows 
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consists of finding functions (■0, h) and the above coupling variables such that 



< 



solves (5) in X [0,T], 
h solves (12) onX X [0,T], 

{ip,h)eA on2:x[o,r], 



(13) 



where A denotes the set of physically admissible states {ip,h}. The admissible 
set A (see Figure 2) has two branches, the branch {h = 0} is associated 
with the dry surface where the soil hydraulic head is less than or equal to 
zero corresponding to unsaturated conditions, while the branch {h = ip} is 
associated with the wet surface where the soil is saturated and the hydraulic 
head is in hydrostatic equilibrium with the overland flow pressure. Thus, the 
admissible set A is defined as 



We mainly focus here on situations where the overland flow is produced by 
exfiltration. Indeed in this situation, a smooth behavior on the admissible set 
can be expected. More drastic situations like runon surface waves on unsat- 
urated soils can in many cases lead to a departure from the admissible set 
especially if the soil is too dry. In these limit situations, other models can be 
more suitable: for instance when inflltration processes are very slow, a model 
where surface flow coexists with an unsaturated soil can be envisaged. 



A= {{ij,h)eR\ h = ip+}, 



(14) 



where = ^{ip + lipl) is the positive part of ip. 



h 



dry 




Figure 2. The admissible set A. 
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3 Discretization 



3. 1 Discretization of Richards ' equation 

Let {7t,}(,>o be a shape- regular family of unstructured meshes of VL consisting 
for simplicity of affine triangles. The meshes can possess hanging nodes. For 
an element r e Tf,, let 9r denote its boundary and its outward unit normal. 
The discontinuous finite element space Vf, is defined as 

{ct> e L'^ip), Vr e r^, 01, e P^lr)}, (15) 

where Pp(t) is the set of polynomials of degree less than or equal to p on an 
element r. Wc observe that the functions in Vf, need not be continuous. This 
fact is exploited by selecting basis functions which are locally supported in a 
single mesh element. The set of mesh faces is partitioned into T\\JT^'^\JJ-'^ 
where ^ is the set of internal faces, JF^^ the set of faces located on W U B 
and the set of faces located on X. For a face F e J^, there are r+ and t~ 
in 7f| such that F = dr^ fl (9r~ and we define the average operator {\p and 
the jump operator as follows: for a function ^ which is possibly two- valued 
on F, 

{e}F = ^(r+n and m^'^^--^, 

where = ^|i-±. For vector- valued functions, average and jump operators 
are defined componentwise. We define to be the unit normal vector to F 
pointing from r~ to . The arbitrariness in the sign of the jump is irrelevant 
in the sequel. 

In the present work, faces on I can exclusively be fiagged either as dry or 
as wet, that is, we do not track the wet/dry interface inside such faces. As a 
result, the set JF^ can be further divided into JF^'* and where JF^'* collects 
the faces flagged as dry and those flagged as wet. These two sets of faces 
induce a partition of T as X[,^'* U X^'*, where 

xj'* {x el,3Fe J^'*, xeF} and T^'* {x el,3Fe J^^'*, X e F}. 

Space discretization 

Let be the discrete approximation of ifj. The symmetric interior penalty 
discontinuous Galerkin method for Richards' equation can be concisely written 

as 

VreT^, V0ePp(T), j^dt[e{i^^)](t> + ar{i^^,i^^,(t>) ^K{i^^,(t>), (16) 
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where for {(, ip, 0) e VJ, x V,, x Pp(t), 



Jt J Ot 



briC, 0) = / V • (K(0^z) + KiC, 0). 



(17) 
(18) 



Here, is the numerical flux associated with the hydrauhc head 



VFe^^, ^(^)l^ 



dcf 



if F e J^i 



■w,t 



if F e J^^^'* u JP-J^^, 



and u{(, ip) the numerical flux associated with u =^ —K(ij:)V'ijj, 



-{i^(C)VV^}^ + TiK^d^'MFfiF if F e 

- X(C) VV' + vKsdF^t/jnn 




if F e jr^'*, 



where r] is a positive parameter (to be taken larger than a minimal threshold 
depending on the shape- regularity of Tf,), Kg the hydraulic conductivity at 
saturation and dp the diameter of the face F which is deflned as the largest 
diameter of the triangle(s) of which F is a face. Observe that for a flow in 
a porous medium with variable conductivity (as in variably saturated flows 
because of the dependence of the conductivity on the hydraulic head), the 
penalty coeflicient at a given interface should scale as the harmonic means of 
the normal hydrauhc conductivity on both parts of the interface, see [25,14]. 
Here, the variations of K are sufficiently mild to use simply the hydraulic 
conductivity at saturation. Furthermore, 6^(^,0) collects the parts of the nu- 
merical fluxes on boundary faces which are independent of ip, namely 



I) 



9 



Summing ar{C, i^, <t>) over all mesh elements yields the global form 

- E / i{K{OW<P}m-np + {K{QV^P}M-nF-r^K,d-^^[ 

- E j^{K{()VcP^-nF + K{C,)VtlJcp-nF-riKsd-F^'il^c^y (19) 



The parameter rj must be chosen large enough to ensure that the form Of, is 
coercive, in the sense that there is a > such that for all G Vf,, 

a^(0») > «( E / i^(0)|V(/>|^+ E / [0f + E ^^^f' / 0')- 



Time discretization 

Let Nt be the total number of time steps and let 5t be the time step taken 
to be constant for the sake of simplicity and such that Nt == T/6t is an 
integer. For any function of time x and for any integer n > 0, denotes the 
value taken by x time n8t. Furthermore, the time derivative of x can be 
approximated by a backward differentiation formula [27] in the form 

i^tXT ^ E 77X"-^ (20) 

r=0 

where q is the order of the formula and {a^}o<r<g are suitable coefficients. 
Using the approximation (20) in (16) for each n e {1 • • ■ Nt} leads to 

Vrer^, V0eP,(T), 

^ I e{r^)<j^ + arir^, r^, 0) = wir^, 0) - E ^ / o{rr)<i^. (21) 

For the first few time steps, a BDF of lower order or a one-step implicit 
scheme can be used, for example a diagonally implicit Runge-Kutta scheme 
or the Crank-Nicolson scheme. The former can present the drawback that the 
last stage can lead to difficulties where the soil is being saturated. 



Nonlinear iterative solver 

The nonhnear equation (21) is solved by the iterative algorithm outhned in 
Algorithm 1. The discrete functions {^^~'^}^^ < being known, successive ap- 
proximations ip^""^ of iph are computed using a quasi- Newton procedure of the 
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form 



Vr e T^, V0 G Pp(r), 

f / (^(V'r) + dAoir,nwr^' - rr))<f> (22) 

Let = tp^'"^'^^ —ip^'"^ and let dj- be defined as dr{C: 4': 0) =^ j 9^[(^{C)]'ip(p, 

so that equation (22) can be written as 

Algorithm 1 Iterative algorithm at each time step for solving Richards' 
equation 

T I n—l in— 2 i^—q / n,0 

Input: ,-0^ ^-^t, ,e„igi 

set m = 
repeat 

Find 5^)1'"" e ^ solving (23) 

set ■^[j — si^f) fA'f) 
m <— m + 1 

until < Ealgl 

Output: ^p^ = V^'"" 



The simplest initialization of Algorithm 1 consisting of choosing the approxi- 
mation of the solution at the previous time step {ip^''^ — V'^"^)) but a higher 
order initialization can also be used (see §4). The error measure E is the rela- 
tive Euclidean norm of the component vector associated with Sip^'™', and Caigi 
is a user-defined convergence criterion. 

3.2 Discretization of the kinematic wave equation 

The kinematic wave equation is discretized on a surface mesh on I which is 
simply the trace of the mesh Tf, on X. Let Nj be the number of mesh faces 
covering I. A finite volume scheme with Godunov flux and time step 6t' is 
used. The time step is taken less than or equal to the time step for Richards' 
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equation, that is, 5t' = St/n' with n' > 1 (see Figure 3). This choice is made 
because the exphcit FV scheme is, as usual, restricted by a CFL condition to 
ensure its stabihty. This is not the case for the discrete Richards' equation 
where, owing to the use of a BDF, a larger time step can be employed. This 
leads to the following notation: h^''^ for n e {1 ■ ■ ■ Nt} and k e {0---n'} 
denotes the discrete approximation of h at time nSt + kSt' and for brevity we 
write = = h^~^'^ . Let Xi, li, x^_i and x^^i be defined on a generic 



2 "^2 



t . , t + n'5t' 

overland flow 



subsurface flow ^i - 



time 



t t + 5t 

Figure 3. Multiple time-stepping for subsurface and overland flows. 

mesh face on I respectively as the center, the length, and the left and right 
vertices of Cj (see Figure 4). Si denotes the slope of the face Cj. Since the flux 
function (/? is convex and the surface water depth is nonnegative, the Godunov 
flux coincides with the upwind flux, yielding 



VA; e {l---n'}, Vi e {1---Nj}, 

1. 



where for all i e {1 • • • Nj}, h^'^ ^ K)''^\ei and v^'" is a discrete interface flux 
yet to be deflned (see §3.3). Observe that a flxed interface flux is used for 
the multiple time steps comprised in a single time step of Richards' equa- 
tion. Equation (24) requires the knowledge of the surface water depth at 




t = (initial condition) and to the left of the first face on a fictitious cell 
(boundary condition) at all discrete times, Vi G {!• ■ -Nj}, = h^{xi) and 
Vn e {1 • • • NT}yk e {0 • • • n' - 1}, Will = h""/. 
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The CFL condition for the exphcit scheme (24) can be expressed as 



where /imax is an a priori bound for the surface surface water depth h on 
X X [0, T]. By definition of the flux function ip, this yields 



In the absence of rainfall and coupling terms, the satisfaction of the CFL 
condition implies a discrete maximum principle and a decrease in the total 
variation for the discrete surface water depth. 



3.3 Single-step coupling algorithm 

We consider in this section the case where Richards' equation is discretized 
in time using a first-order BDF (that is, the Euler implicit scheme). Together 
with the finite volume scheme described in §3.2 for the kinematic wave equa- 
tion, this yields a scheme to approximate the coupled system (13) provided 
we specify the time evolution of the coupling variables {l^'",X^'",ci;",cc;^} for 
n e {1 • • -Nt} (see §2.4). Here, as before, the superscript n stands for the 
value at nSt, so that X^'" = X^'""^^ and so on. This time evolution is designed 
with the twofold objective to ensure that a suitable approximation of {ip, h) 
lies in the admissible set A at all discrete times and to ensure overall mass 
conservation for the whole system (subsurface and overland flows). The re- 
sulting algorithm is outlined in Algorithm 2. It is termed single-step coupling 
algorithm in reference to the use of the first-order BDF which spans a single 
time step interval. For simplicity in the presentation of Algorithm 2, we define 

• V'^ ^ Richards_BDFl(Xj'",X^'",u;^,a;;^, V^-^) as the resolution by Algo- 
rithm 1 of Richards' equation on a time step by the SIPG method, the first- 
order BDF and boundary data on X determined from {X^'",X^'", a;", o;^}, 

• <— Kinematic_wave(/i^~^, n', Ur, v^'") as the resolution of the kinematic 
wave equation by using (24) n' times, 

• v^'^ <— Normal_Velocity(X^'",X^'", cj", cj^, ■0^) as the evaluation of the in- 
terface normal velocity on X defined as 




(25) 




def 



if F e It' 



V 1^ 



(26) 



^v{7jj^\F)-nn + r]KsdF\^^\F 



o;;^!^) ifFeX^'". 
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Note that the expression for v^'" on 2"^'" corresponds to the normal compo- 
nent of the H{dw, 0)-conforming velocity reconstruction derived in [13] for 
DG methods. 



Algorithm 2 Single-step coupling algorithm 

Input: ^r^ 

<— Kinematic jwave{h^~^,n', v^~^, 0) 
Set p = and /i^'^ = ^, 

repeat 

p p+ 1 

Partition of I: Xj'"'^ = {a G J^^,\/k G {1 • • -p}, hJl'^ < 0} and X^'"'^ = 

uy^-^/5t onXj''*'^ 
onX™ 

V'^'^ ^ Richards_BDFl(jJ'"'^ J™,o;r,o;7,^^-i) 
^ Normal_Velocity(jJ-'^ J-"'^<'^a;;'^^^'^ 

VzG{l...iVx},/.r'^ = /.r + 5tA4C'' 
until V2 G ■iVi},/i"'P > 

Output: - V'^'^ and /i^ = /i^'^ 



The principle of Algorithm 2 is the following. Firstly, the surface water depth is 
predicted without subsurface coupling term {v^"' = 0). This predicted surface 
water depth then serves as a Dirichlet boundary condition for Richards' 
equation. Because the Godunov scheme satisfies a discrete maximum principle, 
> for all i G {1 ■ • ■ Aj}, so that jj'"'^ = and J^'"'^ = J. That is, we 
begin the iterations by assuming that X is totally wet. Thus, for p = 1, the 
determination of cu"'^ is irrelevant. Then, Richards' equation is solved and a 
first estimate of the normal velocity v^'"'^ is used to evaluate the surface water 
depth /i^'^. The sign of /i^'^ is subsequently checked on the faces of J. If h^'^ 
is nonnegative on all faces, the evaluation of the hydraulic head and of the 
surface water depth can be accepted as the solution to the coupled system 
at time n5t. Otherwise, a new partition of X is determined and a Neumann 
condition is enforced on those faces where the surface water depth is negative. 
This Neumann condition is evaluated in such a way that at the corresponding 
interface cells, the surface water is completely infiltrated into the soil since 
^v'^ = —h'^/St. A new hydraulic head and a new surface water depth are then 
calculated and the loop is repeated until convergence. Note that convergence 
occurs since the set if '"'^ increases with p while the set xr''"'^ decreases. 
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Admissibility of{ip,h) 



An important point is that Algorithm 2 dehvers nonnegative surface water 
depths. Moreover, on the wet part of the interface, there holds 

since the value of the Dirichlet data o;^'^ on J^'"'^ is fixed during the loop. 
This is not the condition i/j ~ h enforced by the admissible set A but an 0{St) 
approximation of it. Furthermore, on the dry part of the interface, the surface 
water depth is equal to zero and there holds 

yne{i---NT}, VFejJ'", ^;\F<h;\F. 

Again, this is an 0{6t) approximation of the condition < enforced by 
the admissible set. Furthermore, we observe that, if on a given face e^, the 
surface water depth is zero as well as the upwind fluxes over the time 
step [{n — l)St, nM], the Neumann condition on Richards' equation is equal to 
the rainfall intensity. Moreover, we observe that in contrast to front tracking 
schemes. Algorithm 2 does not use any information from the previous time 
step to determine the wet portion of the interface. This offers the advantage 
of robustness and ease of extension to 3D/2D settings, but can entail higher 
computational costs than those incurred by front tracking schemes in the 
absence of exfiltration (see for instance [7]). 

Overall mass conservation 

The total volume of water in the domain Q at time n6t is obtained by inte- 
grating the volumetric water content in Q 

Taking the test function (f) equal to 1 in the SIPG scheme (23), summing 
over the mesh elements and using the first-order BDF to approximate the non 
stationary term yields 

K^nd - V^-^ = {f^ + F^^) 5t + e\ (27) 

where F^ (resp. Fy^^) is the flux over the time step [(n — l)5t, n5t] across the 
interface X (resp. across the bottom and lateral walls), 

F^'^-fvl'- and F^^'M-f v^, (28) 

and e" represents the numerical error in the resolution of the nonlinear sys- 
tem. Recall that |e"| < Ctaigi, where Caigi is the convergence tolerance of 
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Algorithm 1 and C a constant due to the fact that the convergence criterion 
in Algorithm 1 limits the norm of the variation of the hydraulic head ip rather 
than the one of the volumetric water content 9{il)). The total volume of water 
in the overland flow at time n5t is obtained by integrating the surface water 
depth over X 

^over — Jj^' 

The total variation of water volume in the kinematic wave equation over the 
time step [{n — n5t] is obtained by summing the elementary contributions 
in equation (24), yielding 

K'^er - V:-J ={-F^ + Fl^) 6t, (29) 

where Fj is already defined above and where F^^^ represents the water flux 
over the time step [{n — l)St, n5t\ due to the rain and the discharge at points 
A and B, with F^b^ = + F§ + F^, and 

rj./ n' fj./ n' stj./ n' „ 

The total volume of water contained in the coupled system is the sum of the 
volume of each system, V = V^^^^ + V^^^,. When (27) and (29) are summed, 
the interface flux cancels, yielding 

- = (F-e + Fl^;)St + e" (30) 

This relation readily implies the following overall water volume conservation 
result for the single-step algorithm. 

Property 1 Let SV^ be the overall water volume defect over the time step 

[{n - l)5t,nSt] defined as SV" = - 1/"-^ - (F;^^ + F^^^)St. Let AV" 
be the overall water volume defect over the time interval [0, nSt] defined as 

AV'' = J:7=iSV\ Then, 

|A\/"| <nCe,;,i, (31) 
where C is a constant and Caigi is the tolerance in Algorithm L 



3.4 Two-step coupling algorithm 



We consider in this section the case where Richards' equation is discretized in 
time using a second-order BDF for which 

(|)".J_(!,,.-2.-.lv-). (32) 
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The single-step coupling algorithm is not conservative when the non-stationary 
term of Richards' equation is approximated by the second-order BDF owing 
to the fact that the Euler explicit scheme used to solve the kinematic wave 
equation spans only a single time step. Consequently, to obtain a mass con- 
servative scheme, the interface flux used in the kinematic wave equation 
needs to be transformed into a new interface flux so that (29) becomes 

Cer - Ce"/ + ^Br) ^t. (33) 

To identify the expression for observe that using a second-order BDF 
modifies (27) into 



which can be rewritten as 

: (^grnd - V^-^ - \ {y^ll - V^lt) - FwB^t = F^St + 6", (34) 



3 

5' 



where the fluxes F^ and Fy^^ are still defined by (28). Moreover it results from 
(33) that 



2 (K)ver ~ Kwer"*^) ~ 2 (^"^''^ ~ Kver^) ~'~ ( ~ 2^^^^ ~^ 2^^^^^^^ 

-{-l^x + l^r'yt. (35) 



The new interface flux $5 is determined so that the mass flux Fj is exactly 
counter-balanced by the interface flux in (35), whence 



2^2^ J3X3 



At the first time step where a one-step implicit scheme is used, water volume 
conservation is directly enforced by setting =^ Fj. 

The resulting algorithm, referred to as two-step coupling algorithm, is outlined 
in Algorithm 3. Only the differences with Algorithm 2 are indicated. The key 
modification concerns the evaluation of the interface normal velocity in the 
calculation of h^'^. The Neumann data cu"'^ is also modified to ensure that 
the Neumann condition indeed leads to a dry state in the corresponding cell. 
Also, the discrete approximation at time (n — 2)St is added to the input 
and the interface normal velocity v^'"' is added to the output at each time step 
since it is used in the subsequent time step. 

The main result concerning the overall water volume conservation for Algo- 
rithm 3 is the following. 
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Algorithm 3 Two-step coupling algorithm 

Input: ^r^ 

repeat 

^ -(3/i^/5t + t;;'"-')/2 on Jj'"'^' 

vi e {1 • • • TVx}, /^r" = + 5tlk UH''''' + 

until V« e {l---Nx},hl^ > 
Output: ?''|"./'[| and (y" = (y"''^' 

Property 2 Le^ 6e t/ie overall water volume defect over the time step 
[{n — l)St,nSt] defined as 

^yn djf yn _ yn-1 _ ^pn^^ ^ F^^^)5t, (36) 

where Fy^^ =^ ^^wb + ^-^vvb ■ AV"' be the overall water volume defect 
over the time interval [0, nSt] defined as before. Then 

\AV^\ < ^\5V'\+nCeaigi, 

where C is a constant and Caigi is the tolerance in Algorithm 1. 

Proof: Owing to (35), the coupling terms are eliminated when (33) et (34) are 
summed, leading to 

\ {V^ - V^-') - \ iy--^ - V--") - [P^s^t + \fI^^H - ^FXi,^) St^e\ 

Using the definition oi6V"' yields the recurrence relation SV"" = ^SV"'~^ + ^ e"', 
so that 

1 2 " 

3- 3^ 
Owing to the triangle inequality, it is inferred that 

' ' ~ 2 V3 3"+!^' ' ^ V 3'+V' ' 
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whence 

1 n 1 

|Ay"| < -\5V^\ + ^ |e*| < -\5V'\ + nCe„^,i. □ 

Finally, we observe that similar developments can be considered for arbitrary 
order BDFs with additional technicalities and longer recursion formulas. 



4 Results 



Algorithm 3 is assessed on three test cases: the first one concerns overland 
flow over a variable topography, the second one infiltration due to rainfall and 
the third one exfiltration resulting from injected water at the bottom of the 
aquifer. The soil consists of sand and is parameterized by the Haverkamp's 
constitutive relations [17] 



with parameters 



6 s — Or 

1 + |q;^|/ 



9r and K{7P) 



Ks 



i+i^v^r 



= 0.5, a = 0.028cm-^ Ks = IQ-'^cm.s-^, 7 = 4, 
Or = 0.05, /3 = 4, i = 0.030cm-^ 

Figure 5 presents the volumetric water content and the hydraulic conductivity 
as a function of the hydrauhc head. Furthermore, the Strickler coefficient K, 
is set to 60m^/^s~^. For all test cases, the bottom boundary B is located at 
z = 0. 



o 



Hydraulic head ■0 




Hydraulic head ^ 



Figure 5. Hydraulic properties of the soil used in the test cases. 

Piecewise affine finite elements are used {p = 1 in (15)) along with the usual 
local Lagrangian basis functions. For the first time step, the Crank-Nicolson 
scheme is used. A direct solver based on the LU decomposition is employed 
to solve the linear systems. The convergence tolerance eaigi in Algorithm 1 
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is set to 10~^ and the parameter rj is set to 10. Moreover, we focus on the 
use of the second- order couphng algorithm. This choice is motivated by the 
fact that it yields second-order discretization errors in time along with second- 
order discretization errors in space in the L^-norm since p = 1. In addition, a 
second-order initialization of Algorithm 1 is chosen in the form 

Vn>3, i^;^' ^ 3^^-' - 3^^-' + (37) 

2 

except for the second time step where the first-order initialization i/j^^' — 
2'0^ — is used. The second-order initiahzation (37) can decrease significantly 
the CPU time in comparison with the initialization ijj^'^ — 



4.1 Test case 1 (TCI) 



In this first test case, the runoff ffow and the drainage of the subsurface domain 
is induced by the presence of the outlet, located below the initial height of 
the water table. The geometry is presented in Figure 6. The interface X is 
divided into three parts, Xi — {{x,z) e X, x e [0,1.4]} (slope Ji = 0.1%), 
J2 = {{x,z) e J,,T G [1.4,1.6]} (slope J2 = 0.3%), and J3 = {{x,z) G J,a; e 
[1.6,3]} (slope J3 = 0.1%). The final simulation time is T = 300s. The initial 
condition is an horizontal water table located at 0.3025m with an hydrostatic 
pressure profile and the boundary condition on walls and bottom is a zero 
flux, 

= + 0.3025m in O, 

vjv = on (>VU,B) X [0,r]. 

For the overland flow, the initial condition is a horizontal free surface and the 
boundary condition is a zero water depth 

/iO = -z + 0.3025m on /, 

/lA = at A X [0,T]. 



0.3025m 



Jl = 0.1% J2 = 0.3% J3 = 0.1% 

JZ_ 



0.3m 



^ = 



3m 



Figure 6. TCI - Geometry, initial water table position and hx- 
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A mesh with 2063 triangles (corresponding to a typical mesh-size of 3.5cm) 
along with time steps 5t = 2.5s and 6t' = 0.25s have been used. We have 
verified that the interface normal velocity obtained with 6t = 6t' = 0.25s can 
be superimposed to that reported below. In this case, the use of St — 2.5s 
instead of St — St' — 0.25s leads to a gain of 89% in the computational 
time (i.e. the computation effort required for performing one time-step in the 
solution of Richards' equation is eighty times more expensive than the one of 
the kinematic wave). 

Figure 7 presents the free surface of the overland flow (/i^ + topography) and 
the interface normal velocity v^'^ along the interface at three characteristic 
times of the simulation (10s, 100s and 300s). The free surface being piecewise 
constant, it is depicted on each interface cell by a solid line. The interface 
normal velocity v^'"' is plotted with circles if the interface is wet (that is, on 

J^' '") and with crosses if the interface is dry (that is, on X^'"). 

Figure 8 provides a closer insight at the issue of staying on the admissible set 
A. For the same times as in Figure 7 and for each face of JF^, each couple 
[ipl^, h^) is represented by a cross (the mean-value of ipl^ is considered on each 
face). The admissible set A is also plotted with a solid line. 

The hydraulic jump in the overland flow is visible at the beginning of the 
simulation on Figure 7 at 10s. Moreover, exfiltration appears on some faces 
located on X2 and X3. During the simulation, a Neumann boundary condition is 
imposed on the faces where the water becomes equal to zero. It is confirmed by 
Figure 8 where the number of couple h''^) situated on the branch {h = 0} 
increases. 



4.2 Test case 2 (TC2) 



The principle of this test case is inspired by the work of Abdul and Gilham [1]: 
a constant rainfall intensity is imposed on the upper part of the domain for a 
fixed period of time, whereas the lateral and lower boundaries are imperme- 
able. In our case, the geometry is shown in Figure 9 and the final simulation 
time is T = 360s. The initial condition is an horizontal water table located 
at 0.85m with an hydrostatic pressure profile and the boundary condition on 
walls and bottom is a zero flux, 

tlP ^-z + 0.85m in Q, 

vn^^ on (WUB) X [0,r]. 
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Figure 7. TCI - Free surface (solid 
line) and interface normal velocity 
{cm/12min) plotted with circles if in- 
terface is wet and with crosses if inter- 
face is dry. 




0.1 0.2 



Figure 8. TCI - Cloud of points 
(V'^, h"^) on the admissible set A at dif- 
ferent times. 



For the overland flow, the initial condition and the boundary condition are 

/i° =0 on I, 

/lA = at A X [0,T]. 

A constant rainfall intensity equal to 10% of the hydraulic conductivity at 
saturation is imposed during 180s and is stopped afterwards, 



Vr-m^-OAKs onXx [0,180], 
v,-nn^O onXx [180,r]. 
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O.IK^ = 3.6mmh ^ 



0.85m 




6m 



Figure 9. TC2 - Geometry, initial water table position and constant rainfall intensity. 

A mesh with 2049 triangles (corresponding to a typical mesh-size of 10cm) and 
time step 6t = 5t' = Is have been used. We have verified that the interface 
normal velocity obtained with a finer mesh (8763 elements) and a smaller 
time step {St — 5t' — 0.5s) can be superimposed to that reported below. Also, 
observe that 5t' — Is roughly corresponds to the CFL condition, so that, 
for the present test case, the accuracy limit on the time step for Richards' 
equation is comparable to the CFL restriction. 





100 200 300 400 500 600 100 200 300 400 500 600 



Figure 10. TC2 - Free surface (solid line) and interface normal velocity {cm/lbmin) 
plotted with circles if interface is wet and with crosses if interface is dry. 

Figure 10 presents the free surface and the normal velocity f along the in- 
terface at four characteristic times of the simulation (10s, 60s, 180s and 360s). 
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The same notation is used as in Figure 7. The hydrological response of the 

system can be divided into four phases. 

1 - Soil saturation [0, 50s]. In this phase, which results from the initial water 
table position, the 15cm top layer is being saturated. The rainfall is totally 
absorbed by the soil, the surface water depth is equal to zero and a Neumann 
condition is imposed on the all faces of 1 (Figure 10 at 10s). 

2 - Surface runoff occurs on part of I [50s, 90s]. The rainfall is still partially 
absorbed by the soil, but a Dirichlet condition is now being imposed on the 
part of X located near the outlet. Interestingly, infiltration occurs on the most 
part of the interface since the normal velocity is negative, but some exfiltration 
occurs on the first few faces located near the outlet where the normal velocity 
becomes positive (Figure 10 at 60s). 

3 - Surface runoff occurs on X [90s, 180s]. Surface runoff occurs on the whole 
interface and the soil is totally saturated. A Dirichlet condition is imposed 
throughout the interface and the surface water depth is positive (Figure 10 at 
180s). 

4 - Drainage [180s, 360s] . When rainfall stops, the surface water depth returns 
to zero on the faces located near the point A because of infiltration and surface 
runoff. A Neumann boundary condition is imposed on the dry zone near the 
point A (Figure 10 at 360s). 

Figure 11 provides a closer insight at the issue of staying on the admissible 
set A. For the same times as in Figure 10 and for each face of each couple 
(-0^, h^) is represented as in Figure 8. Note that different scales are used, 
so that the branch {h — ■0} is almost vertical. The four phases described 
above are clearly illustrated by the position of the cloud of points. At 10s, the 
hydraulic head is negative and the water depth is equal to zero. The cloud of 
points is only on the branch {h = 0} corresponding to a dry soil. At 60s, the 
hydraulic head is equal to the water depth for some faces. The cloud of points 
is located on the two branches because the soil contains both saturated and 
unsaturated zones. At 180s, the hydraulic head is equal to the water depth for 
all the faces. The cloud of points is only on the branch {h = ip} corresponding 
to a wet soil. At 360s, the hydraulic head becomes again negative where the 
surface water depth is equal to zero. The cloud of points is again located on 
the two branches. 



Figure 12 presents results on mass conservation issues. Multiplying equation 
(36) by the water density p, summing over the time intervals in [0, n6t], know- 
ing that -F>vgn = and the definition of F", F^^^ and AF" yields 



AM" is defined as the total mass variation over the time interval [0,n(5t] 
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Figure 11. TC2 - Cloud of points (V'^, fi^) on the admissible set A at different times. 

and is the sum of the total groundwater mass variation AM" and the to- 



tal overland mass variation AM"^gj,. The quantities Th=i EIL 



out 



and 



are respectively the total inflow of water, the total outflow of water and 
the total mass balance defect cumulated at time nSt. The five quantities 
AM", AM^,„d>^M;jver>Er=iM^„ and EILi M^Jut are presented in the left part 
of Figure 12. In particular, this figure confirms the four phases of the sim- 
ulation. The rainfall is totally absorbed by the soil at the beginning of the 
simulation until 50s since AM" = AM^^^^. Then, the increase of AM^^^ di- 
minishes and AM"ypj. becomes positive as a result of soil saturation. From 90s 
to the end of the simulation, the variations of AM" and AM^^^jj, are the same, 
corresponding to a complete saturation of the soil. Moreover, during the last 
phase, the total water inflow is constant because the rainfall stops, so that the 
total water outflow is the same as the total mass variation. 



Both total mass balance defects obtained with the single-step and two-step 
coupling algorithms are compared in the right part of Figure 12. While Algo- 
rithm 3 yields a sizable improvement over Algorithm 2, it can still be noticed 
that the mass balance defect produced by Algorithm 2 is only of the order of 
a few percent of the global quantities such as AM". 

Finally, Figure 13 studies in more detail the mass fluxes in the kinematic 
wave equation. The mass flux is decomposed into the exflltration flux 
and the inflltration flux in the form = + with 

pn,+ <M _ j^^^ ^*,n pn- d£ _ _ ^*,n^ ^-^^^^ time-dependent sets 
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Figure 12. TC2 - Left: Mass repartition in the coupled system; Right: Mass balance 
defect E" for Algorithm 2 (dashed) and Algorithm 3 (solid). 

2^''^ and 2^'" are defined as follows 

J^'+ {x e I; vl'^'ix) < 0} and J^'' {x e I; vl'^'ix) > 0}. 

The four quantities pStFj'^ , p5tF^~ ,Mf^ and M^^^ are plotted on Figure 13 
as a function of time. 
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Figure 13. TC2 - Mass fluxes in the kinematic wave equation. 



4.3 Test case 3 (TC3) 



In this third test case, an exfiltration is produced on the upper part of the 
domain by injecting water at the bottom-left part. The geometry is presented 
in Figure 14 and the final simulation time is T = 360s. The initial condition is 
an horizontal water table located at 0.1m with an hydrostatic pressure profile 
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and the boundary condition on the walls corresponds to a zero flux, 



i/j'^ ^ -z + O.lm in f], 

vn^O onyVx[0,T]. 

The rainfall intensity is set to zero. An infiltration fiux with a parabolic profile 
and a mean- value vn equal to 5% of hydraulic conductivity at saturation is 
imposed during 2 minutes on the left half Bi of the bottom [Bi = {{x,z) G 
B,x e [0, 1]} and Br — {{x,z) e B,x e [1, 2]}). This infiltration flux is linear 
during the first 10s, constant on [10, 120], and equal to zero for t > 120s: 



VN{x,t) 



x{x - 1) 0.003^5 t, if {x, t)eBiX [0, 10], 
x{x - 1) OmK,, if {x, t)eBiX [10, 120] , 
0, otherwise. 



For the overland flow, the initial condition and the boundary condition are 

/i° =0 onX, 

/lA = at A X [0,T]. 

A mesh with 1874 triangles (corresponding to a typical mesh-size of 2.5cm) 
and time step 6t = 6t' = Is have been used. We have verified that the interface 
normal velocity obtained with a finer mesh (7310 elements) and a smaller time 
step {5t — 5t' — 0.5s) can be superimposed to that reported below. 

J = 0.2% 



0.1m 



Figure 14. TC3 
groundwater. 
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Geometry, initial water table position and flux infiltration in 



Figure 15 presents the free surface and the normal velocity f^'" along the 
interface at six characteristic times of the simulation (5s, 35s, 50s, 100s, 150s 
and 360s) and Figure 16 presents the surface water depth at these different 
times. The same notation is used as in Figure 7. The hydrological response of 
the system can be divided into six phases. 

1 - Soil saturation [0,15s]. This phase results from the initial water table 
position. At the beginning of the simulation, the soil is partially saturated and 
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Figure 16. TC3 - Surface water depth at different times. 
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the injection at the bottom of the domain increases the hydrauhc head. The 
interface is totally dry and a zero Neumann boundary condition is enforced 
everywhere on X (Figure 15 at 5s). 

2 - Partial exfiltration [15s, 45s]. The soil becomes saturated in the left part of 
the domain and the interface normal velocity positive, so that water begins to 
exfiltrate from the faces situated in this saturated zone. A Dirichlet condition 
is enforced on those faces (Figure 15 at 35s). 

3 - Full exfiltration [45, 100s]. When the soil is totally saturated, the amount 
of exfiltrated water is equal to the amount of injected water. We observe that 
overland flow occurs over the whole interface X and that a Dirichlet condition 
is being enforced everywhere. However, most of the overland flow still remains 
concentrated near the upper part of the interface (Figure 15 at 50s). 

^ - Propagation of the runon wave [100, 120s]. In this phase, the runon wave 
propagates downstream. It is worthwhile to notice that a slight part of the 
surface water infiltrates back into the soil as indicated by the sign of the 
normal velocity near the heading part of the runon wave (Figure 15 at 100s). 

5 - Outflow [120,200s]. When water injection ceases at the bottom of the 
domain, the amount of exfiltrated water decreases sharply and there is even a 
small portion of the interface located near the point A where water infiltrates 
back into the soil (despite the boundary condition is of Dirichlet type since the 
surface water depth is still positive) while most of the overland flow reaches 
the outlet and exits the system (Figure 15 at 150s). 

6 - Drainage [200,360s]. The surface water depth vanishes on the upper part 
of the interface X where a zero Neumann condition is now imposed (Figure 15 
at 360s). 

Figure 17 shows that each couple (z/'^, h^) stays on the admissible set A. As 
in the previous test case, the phases described above are clearly illustrated by 
the position of the cloud of points. It is located on the branch {h — 0} when 
the soil is unsatured, on the branch {h — ■0} when the soil is saturated and 
on the two branches when there are both saturated and unsaturated zones at 
the interface. 

Results on Figure 18 and Figure 19 are similar to the ones of the previous test 
case, in particular the comparison of the mass balance defects for the one-step 
and the two-step algorithms. 



5 Conclusion 

In this work we have presented a robust and accurate numerical method to sim- 
ulate coupled subsurface and overland flows governed by Richards' equation 
and the kinematic wave equation. Special care was taken to design coupling 
algorithms that preserve the overall mass in the system and that also satisfy 
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the various equality and inequality constraints imposed at the interface. Ex- 
tensions of this work include the use of more complex models, such as the 
shallow-water equations, to describe the overland flow and the possibility to 
account for drainage pipes in the soil. Extension to two-dimensional surface 
flows and three-dimensional subsurface variably saturated flows can also be 
considered. The present algorithms are currently being tested in more complex 
and realistic test cases related to held studies. 
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Figure 18. TC3 - Left: Mass repartition in the coupled system; Right: Mass balance 
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MP 


4 




in 


2 






ass 


\ \ 






-2 




p5tF~^^^ e ^ 


-4 




^'rout , , , 



Time 

Figure 19. TC3 - Mass fluxes in the kinematic wave equation. 
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